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ABSTRACT 

Shape  structural  optimization  for  blast  mitigation  seeks  to  counteract  the  damaging  effects  of  an 
impulsive  threat  on  occupants  and  critical  components.  The  purpose  of  a  vehicle  energy-deflecting 
hull  is  to  mitigate  blast  energies  by  channeling  blast  products  and  high-pressure  fluids  away  from 
the  target  structure.  Designs  of  pyramid-shaped  protective  structures  have  been  proposed  in 
existing  and  concept  vehicle's  platforms.  These  structures  are  more  effective  than  the  traditional 
flat-plate  designs  in  terms  of  cabin  penetration  and  weight.  Studies  on  other  blast  concept  design 
remain  scarce.  This  investigation  addresses  the  design  of  blast-protective  structures  from  the 
design  optimization  perspective.  The  design  problem  is  stated  as  to  finding  the  optimum  shape  of 
the  protective  shell  of  minimum  mass  satisfying  a  deformation  and  envelops  constraints. 
Performance  improvements  are  observed  as  the  envelope  constraint  is  relaxed  and  the 
optimization  problem  includes  a  larger  number  of  design  variables  with  the  consequent 
computational  cost.  In  order  to  consider  the  solution  of  a  problem  with  no  envelope  constraints, 
this  work  explores  heuristic  alternatives:  inverted  profile  and  hybrid  cellular  automata  (HCA). 
Numerical  results  demonstrate  that  convergent  designs  based  on  trigonometric  functions  yielded  a 
greater  reduction  in  mass  over  baseline  than  all  other  candidate  designs,  and  show  promise  for 
future  development.  Designs  generated  via  the  HCA  topography  method  resulted  in  structures 
whose  performance  paralleled  or  exceeded  most  geometrically  constrained  designs  despite  the  fact 
that  they  can  be  created  with  fewer  computational  resources  and  adapted  to  irregular  design 
domains. 
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Trigonometric  function  coefficients 

Element  thickness 

Vector  of  design  variables 

Vector  of  error  values 

Trigonometric  function  coefficients 

Structure  height 

Nodal  location  after  blast 

Radial  coordinate  of  surface  of  revolution 

Time  variable 

Vector  of  x -coordinate  values  for  all  nodes 

Vector  of  y-coordinate  values  for  all  nodes 

Deflection  after  blast 

Vector  of  z-coordinate  values  for  all  nodes 

Polynomial  function  coefficients 

Mass  function 

Number  of  free  nodes  in  the  design  space 
Cabin  penetration  function 
Maximum  allowable  penetration 
Radius  of  surface  of  revolution 
Envelope  constraint  function 
Gaussian  spread  variable 
Gaussian  height  variable 
Element  density 


1  INTRODUCTION 


This  article  investigates  the  automated  design  process  of  isotropic  plates  to  mitigate  the  effects  of 
blast  loading.  Framed  against  the  problem  of  blast  protection  system  design  for  vehicle  applications, 
the  performance  of  the  designs  is  measured  in  terms  of  mass  and  deflection,  as  has  been  done  in 
previous  investigations  of  lightweight  armor  design  [1].  There  has  been  substantial  work  on 
analysis  of  the  effects  of  blast  loading  on  shell  structures.  Argod  and  Belegunda  have  shown 
significant  improvement  of  structure  design  using  velocity-field  based  optimization  and  have 
demonstrated  the  effects  of  different  boundary  conditions  for  plates  of  this  kind  [2].  Methods  of 
blast  energy  absorption  have  been  evaluated  through  extensive  design  investigations  of  composite 
materials  [3]  [4] ,  as  well  as  numerical  [5]  and  experimental  [6]  simulations  of  sacrificial  structures. 

Physical  nonlinearities  are  abundant  in  the  process  of  a  blast  event.  In  such  events,  the  kinetic 
energy  is  transferred  to  the  structure  by  means  of  a  fluid  structure  interaction.  The  dynamics  of  the 
events  and  their  interactions  with  solid  structures  is  an  area  that  has  been  extensively  researched 
[7].  The  governing  equations  of  thermo-mechanics  that  characterize  blast  loading  conditions 
require  the  conservation  of  mass  and  energy,  balance  of  linear  and  angular  momentum,  heat 
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equations,  stress/strain  relations,  and,  in  the  case  of  blast  pressure  calculations  and  shrapnel 
penetration  calculations,  kinetic  and  caloric  equations  of  state;  Kinetic  and  caloric  equations  of 
state  are  required  in  solving  the  behavior  of  pressurized  gases  upon  detonation,  volumetric 
equations  of  state  may  be  used  to  solve  for  the  behavior  of  the  target  ceramic  material.  From  this 
basic  thought  process,  it  appears  that  in  order  to  solve  for  the  behavior  of  the  target  material  while 
accounting  for  all  the  thermo-mechanical  interactions  of  the  blast  wave  and  projectile  impacts,  the 
solution  of  all  sixteen  equations  are  required.  To  reduce  the  total  governing  equations  required,  it 
may  be  necessary  to  simplify  or  dissect  the  phenomena  of  detonation  into  more  manageable 
models. 

Several  experimental  investigations  have  been  conducted  in  the  field  of  detonation  events  [8].  The 
CONWEP  model  developed  by  the  Army  Research  Laboratory,  and  currently  used  in  commercial 
finite  element  blast  simulations,  is  considered  adequate  for  use  in  engineering  studies  of  vehicle 
response  to  the  blast  from  land  mines  [7].  The  CONWEP  algorithm  does  not  explicitly  simulate  the 
effects  of  the  detonation  reaction  with  air  nor  does  it  calculate  the  shock  wave  propagation  and  its 
reactions  with  the  target  structure.  In  continuum  mechanics  the  Arbitrary  Lagrange-Eulerian 
method  is  commonly  applied  in  the  finite  element  analysis  of  fluid  structure  interactions.  The  ALE 
method  is  popular  in  its  application  to  large  shear  deformations  where  a  traditional  Lagrangian 
mesh  breaks  down  due  to  geometric  instabilities.  Compared  to  a  traditional  Eulerian  mesh  such  as 
those  used  in  computational  fluid  dynamics  analysis,  the  ALE  method  is  better  suited  to  tracking 
fluid  boundaries  and  resolving  flow  details  involving  large  volume  changes  [9].  Due  to  the 
importance  of  fluid  structure  interactions  in  many  engineering  applications,  the  ALE  algorithm  has 
been  applied  to  commercial  finite  element  solvers  and  is  currently  available  in  the  LS-DYNA  971 
release. 

The  length  and  time  scales  in  blast  analysis  are  very  different  from  that  of  the  traditional  linear 
elastic  quasi-static  optimization  problems.  In  such  previously  studied  problems  the  material  stress 
strain  relations  are  assumed  to  follow  linear  elastic  models  for  quasi-static  problems  and  linear 
plastic  deformation  relations  are  assumed  for  material  yielding  under  dynamic  simulations. 
Materials  in  the  design  domain  are  generally  assumed  to  be  isometric  and  the  density  of  the 
material  is  assumed  to  be  controllable  and  linearly  proportional  to  the  material  properties.  The 
blast  optimization  problem  proposed  involves  pressure  and  mechanical  loading  conditions  that 
cannot  be  properly  modeled  using  the  existing  linear  elastic  material  models.  The  time  scale 
involved  under  blast  loading  conditions  involves  irreversible  processes  that  require  a  more 
complete  model  of  thermal  mechanical  reactions  than  that  currently  available.  Of  primary  concern 
in  the  detonation  and  explosion  interaction  event  as  opposed  to  a  crash  event  is  the  interaction  of 
the  blast  shock  wave  with  the  target  structure. 

The  fundamental  optimization  problem  to  be  analyzed  in  the  design  of  a  structure  for  blast 
mitigation  is  that  of  minimizing  the  kinetic  energy  transfer  from  the  blast  wave  to  the  solid  body. 
The  goal  of  such  efforts  is  to  develop  a  system  that  could  absorb  a  significant  amount  of  the  energy 
released  in  a  blast  event  such  that  the  underlying  structure  may  be  preserved.  This  article  will 
instead  focus  on  various  design  techniques  for  isotropic  plate  structures  within  the  same  design 


UNCLASSIFIED  Dist  A.  Approved  for  public  Release. 


3 


space,  and  subject  to  the  same  objective  function  and  constraints  in  order  to  compare  design 
methodologies  in  terms  of  optimized  structures  produced  via  formal  optimization  techniques  such 
as  sequential  quadratic  programing  (SQP)  and  hybrid  cellular  automata  (HCA).  The  objective  of  this 
investigation  is  to  compare  feasible  design  methodologies  through  the  expansion  of  the  problem 
dimension  in  order  to  reach  the  limits  of  performance.  For  the  purposes  of  this  study,  seven  profiles 
are  evaluated,  each  under  the  same  loading  and  boundary  conditions. 

The  candidate  profiles  are  divided  into  two  broad  categories:  geometrically  constrained  design 
(GCD)  and  free  shape  design  (FSD).  A  GCD  includes  an  analytical  description  of  the  structure's 
profile  in  the  problem  statement,  while  FSD  such  expression  such  constraint  is  relaxed.  The  small 
number  of  design  variables  describing  the  structure's  profile  in  a  GCD  allows  the  use  of  traditional 
programming  methods  such  as  sequential  quadratic  programming.  In  the  case  of  FSD  the  use  of 
such  traditional  methods  is  impractical;  therefore,  this  work  makes  use  of  heuristic  approaches,  i.e., 
inverted  profile  and  hybrid  cellular  automata. 

2  PROBLEM  DEFINITION 

In  blast  protection  system  design  for  vehicle  applications,  there  are  two  primary  performance 
measures  of  critical  relevance:  weight  and  cabin  penetration.  Weight  is  privileged  as  a  key 
performance  measure  stemming  from  the  need  for  lightweight,  compact  structures  that  can  be 
fitted  to  existing  vehicle  designs  without  extensive  re-design  of  the  frame  or  other  vehicle  systems. 
The  goal  is  to  design  light  structures  that  do  not  adversely  affect  vehicle  performance.  Large 
deflections  of  the  plate  structure  due  to  a  blast  event  can  cause  penetration  into  the  passenger 
cabin  of  the  vehicle  and  result  in  occupant  injury.  In  this  way,  cabin  penetration  can  be  seen  as  a 
means  of  quantifying  the  degree  to  which  the  blast  energy  has  been  mitigated;  if  the  deflection  of 
the  plate  structure  post-blast  exceeds  a  certain  amount,  the  design  is  unsuccessful.  In  consideration 
of  these  performance  measures,  the  objective  function  is  formulated  as  the  mass  of  the  plate 
structure  and  is  minimized  subject  to  displacement  constraints  (cabin  penetration)  and  design 
space  limitations.  The  general  optimization  problem  addressed  in  this  paper  is 

find  d 

minimize  M(d) 

subject  to  Pc(d)  -  Pc  max  ^  0  (1) 

S(x,y,z,d )  =  0 
dL  <  d  <  du 

where  d  E  Rn  is  the  set  of  all  design  variables  characterizing  the  shape  and  thickness  of  the  plate, 
M(d)  is  the  plate's  mass,  Pc(d)  is  the  penetration  after  the  blast  event  with  respect  to  datum  plane, 
Pc  max  is  the  maximum  allowable  value  for  penetration.  The  envelope  constraint  h(x,  y,  z,  d)  is  a 
function  of  the  nodal  coordinates  x,  y,  and  z.  The  box  constraint  dL  and  du  are  the  lower  and  upper 
bounds  for  the  design  variables,  respectively.  The  envelope  constraint  h  is  progressive  relaxed  so 
the  design  space  is  expanded  to  contain  more  design  variables.  This  allows  increasing  the 
performance  design  problem  at  expenses  of  more  complex  topographies.  This  work  is  framed  in  the 
context  of  vehicle  protection  from  blast  events.  Five  envelope  constraints  are  considered  in  this 
work  (Table  1). 
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Table  1:  Envelope  constraints  and  corresponding  number  of  design  variables 
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•  Flat  Plate  Design.  This  topography  for  flat  plane  is  defined  by  the  condition  in  which  the  all  the 
z  —coordinate  values  for  every  node  is  the  same.  This  design  is  regarded  as  the  baseline  for 
comparison  with  each  progressive  candidate  design. 

•  Pyramid  Profile  Design.  Methods  of  blast  mitigating  structure  design  have  been  evaluated 
extensively  through  experimental  efforts  in  the  development  of  lightweight  V-shaped  hulls,  and 
have  demonstrated  that  V-shaped  designs  can  mitigate  the  effects  of  blast  events  [10].  The 
pyramid  profile  design  is  an  improvement  of  the  V-shape  design  used  in  concept  designs,  due  to 
the  fact  that  it  is  constrained  on  four  sides. 

•  Gaussian  Function  Design.  In  the  interest  of  creating  more  complex  curves  through  minimal 
expansion  of  the  design  domain,  a  plate  structure  that  takes  on  the  shape  of  a  Gaussian  function 
of  two  variables  is  generated.  Previous  explorations  of  this  shape  demonstrated  promising,  yet 
inconclusive  results  [11]. 

•  Polynomial  Function  Design.  And  Trigonometric  Function  Design.  Two  additional  design 
methodologies  are  examined  which  relax  the  problem  further.  In  both  the  polynomial  function 
and  trigonometric  function  cases,  a  function  of  several  variables  is  used  to  generate  a  complex 
curve.  The  plate  design  is  achieved  as  a  surface  of  revolution  by  rotating  the  curve  about  the 

z  —axis. 

As  the  number  of  design  variables  increases,  so  does  the  complexity  of  the  analytical  description  of 
each  design.  A  large  increase  in  the  number  of  design  variables  makes  the  problem  intractable  in 
that  an  analytical  solution  for  sensitivity  has  not  been  derived.  In  short,  for  traditional  sensitivity- 
based  design  optimization,  the  cost  of  sensitivity  limits  the  number  of  design  variables.  In  order  to 
overcome  those  limitations,  two  alternative  methods  are  examined  which  utilize  a  large  design 
domain  and  rely  on  user  input  to  avoid  intractability.  These  free  shape  designs  deal  with  the  most 
relaxed  problem,  as  there  is  no  underlying  analytical  description  of  the  shape: 

•  Inverted  Profile  Design.  There  has  been  some  interest  in  creating  a  plate  that  takes  on  an 
inversion  of  the  topography  of  a  flat  plate  after  a  blast  event,  largely  inspired  by  the 
mathematical  idea  of  the  catenary  curve,  the  idealized  shape  that  a  chain  assumes  under  its  own 
weight.  The  reasoning  being  that  the  complex  blast  pressure  distribution  would  then  determine 
the  shape  of  the  plate  and  influence  the  z-coordinate  at  each  node  proportional  to  the 
magnitude  of  the  loading.  This  approach  has  been  utilized  in  the  design  of  bridges  as  well  as  in 
the  world  of  architectural  design  [12, 13]. 

•  HCA  Topography  Design.  The  Hybrid  cellular  automata  method,  or  HCA,  is  an  approach  that 
allows  us  to  handle  many  design  variables.  In  this  approach,  the  field  variable,  in  this  case  the 
nodal  z  —coordinate,  is  driven  toward  a  pre-determined  set  point.  HCA  methods  have  been  used 
in  the  past  by  Goetz  and  Tovar  to  develop  two-material  topologies  for  blast  mitigation  [1] ,  but 
there  is  a  need  to  demonstrate  the  uses  of  HCA  for  topography  optimization. 
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MODELING  CONSIDERATIONS 


2.1  Design  Domain 

The  design  domain,  which  defines  the  envelope  constraints  and  boundary  conditions  for  the 
problem,  is  chosen  is  such  a  way  as  to  simulate  under-vehicle  conditions.  Consider  the  vehicle 
model  shown  in  Figure  1. 


Figure  1:  Truck  model  for  design  space  development  (top)  and  detail  of  under-vehicle  showing  datum 

plane  (bottom) 


The  design  domain  is  chosen  as  a  three  dimensional  space  of  the  same  dimensional  order  as  the 
vehicle  model  shown  in  the  figure,  where  the  datum  plane  (shown  as  the  dotted  line)  is  an  arbitrary 
distance  S  =  0.03  m  from  the  upper  plane  of  the  design  space  to  account  for  the  distance  between 
an  under-vehicle  plate  structure  and  the  passenger  cabin.  This  distance  allows  for  some  penetration 
below  the  datum  plane  without  penetrating  the  cabin. 


In  order  to  privilege  the  display  of  the  optimized  plate  designs,  the  design  domain  has  been 
modeled  from  beneath  or  "upside  down"  orientation,  resulting  in  a  1  m  by  1  m  by  0.15  m  three- 
dimensional  space  with  the  datum  plane  distance  8  from  the  bottom  plane  as  shown  in  Figure  2. 
For  this  orientation,  a  node  is  considered  to  penetrate  the  cabin  if  its  nodal  z  —coordinate  is  less 
than  —0.03  m. 


x 


Figure  2:  Design  domain  representation  with  datum  plane  (z  =  0)  shown. 
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The  penetration  it  is  taken  as  the  z-displacement  relative  to  the  initial  position.  The  thickness  of  the 
plate,  c,  is  computed  as  the  sum  of  the  distance  between  the  through-integration  points  for  the  shell 
elements.  For  all  models,  shell  elements  with  five  through  thickness  integration  points  are  used.  In 
some  design  candidate  cases,  the  nodal  locations  are  determined  as  a  function  of  additional  design 
parameters.  Additional  design  parameters  are  shown  on  a  case-by-case  basis  below. 

2.2  Baseline  Design  and  Mesh  Refinement 

The  baseline  design,  which  serves  a  starting  shape  for  the  optimization  methods,  is  a  flat,  isotropic 
steel  plate  with  the  following  characteristics: 

•  Uniform  thickness:  each  element  has  the  same  through  thickness  based  on  the  distance 
between  integration  points. 

•  Fixed  on  all  sides:  all  degrees  of  freedom  are  constrained  for  all  nodes  at  the  sides  of  the  plate. 
These  boundary  conditions  have  been  demonstrated  to  have  quantitatively  the  same 

z  —displacement  response  over  time  as  a  free  plate  with  a  stiffener  support  around  the  edge  [2]. 
This  is  meant  to  simulate  a  plate  capable  of  maintaining  rigidity  at  the  edges  and  approximate 
conditions  in  vehicle  design. 

The  plate  is  modeled  using  four  node  shell  elements  with  five  through  thickness  integration  points. 
Mesh  density  studies  showed  that  the  z  —displacement  was  relatively  insensitive  to  mesh 
refinement.  In  order  to  preserve  computational  efficiency,  the  plate  was  modeled  with  a  coarse 
mesh  grid  of  26  by  26  nodes,  for  a  total  of  625  elements,  as  shown  in  Figure  3. 


Figure  3:  Baseline  design  26  x  26  node  plate. 


The  plate  is  modeled  as  A36  type  structural  steel  with  linear  elastic-plastic  behavior  -  LS-DYNA 
input  card  *MAT_PIECEWISE_LINEAR_PLACTICITY.  The  following  material  properties  are  used  in 
each  design  case:  mass  density  7800  kg/m3,  Young's  modulus  200  GPa,  Poisson's  ratio  0.3,  yield 
stress  220  MPa,  tangent  modulus  2  GPa. 

2.3  Loading  Conditions  and  Time  Considerations 

The  loading  conditions  were  applied  using  the  ConWep  BLAST_ENHANCED  engineering  model  of  an 
air  blast,  with  the  same  set  of  parameters  used  consistently  throughout  to  ensure  continuity  for 
comparative  purposes:  charge  magnitude  =  5  kg  TNT,  charge  location  =  40  cm  above  the  center  of 
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the  plate  (0,  0,  0.4),  type  of  blast  source  =  spherical  free-air  burst  (LS-DYNA  keyword  input  - 
BLAST=2).  The  peak  pressure  loading  P0,  as  shown  in  Figure  4,  occurs  very  quickly  after  the  initial 
detonation  time  and  the  maximum  z-displacement  occurs  as  a  result  of  this  pressure  [14].  Thus,  a 
very  small  simulation  time  (t  <  0.0055  s)  was  chosen.  Exact  simulation  time  for  each  design  was 
chosen  on  a  case-by-case  basis. 

p 


Positive  Negative 

Phase  Phase 

Figure  4:  Blast  pressure  time  history. 


2.4  Adaptation  of  Optimization  Problem 

Each  candidate  design  is  considered  individually  and  adapted  to  the  existing  optimization  problem 
as  the  number  of  design  variables  in  x  increases.  In  each  case  the  objective  function  and  the  bounds 
on  d  are  different,  but  the  penetration  constraints  and  envelope  are  the  same  case  in  each: 

•  Penetration  constraint:  Pc(d,  t)  —  Pcmax  ^0 

Pc  max  =  0.03  is  chosen  to  ensure  that  the  optimized  topography  results  in  a  design  where  the 
lowest  z  —coordinate  after  the  blast  event  is  -0.03  m,  i.e.,  there  is  no  cabin  penetration. 

•  Design  domain  constraint:  zL  <  z  <  zu 

zL  =  —0.03  m  and  zu  =  0.12  m  are  chosen  to  ensure  that  the  design  does  not  violate  the 
domain  constraints. 

2.5  Approximation  of  Surface  Area 

For  a  given  surface  of  two  variables,  the  area  can  be  determined  analytically  according  to  the 
following  equation: 

ff]1+(£fCx'yrf+(£;fCx'y))  (2) 

For  complicated  surface  geometry  whose  functions  involve  non-linear  equations  of  several 
variables  as  the  ones  examined  here,  analytical  methods  of  determining  the  surface  area  of  the  plate 
become  problematic.  Instead,  the  surface  area  of  the  plate  is  approximated  based  on  nodal  locations. 
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If  we  consider  a  single  element  as  shown  below  in  Figure  5,  we  can  define  vectors  that  begin  and 
end  on  the  nodes  associated  with  each  element. 


node  1  vt  node  2 

Figure  5:  Vector  numbering  to  determine  element  area. 

For  a  given  surface  divided  into  N  elements,  the  surface  area  can  be  described  as  in  (3)  where  a  new 
set  of  vectors  is  generated  for  each  element  (k)  based  on  nodal  coordinate  information. 

A  =  T  \  ||  v\  X  V%\\  +  \  ||  w3k  X  t>4  ||  (3) 

*—lk= l  z  z 

Mass  calculations  based  on  this  method  of  approximating  surface  area  were  compared  to  the  mass 
results  generated  by  the  LS-DYNA  internal  method  and  were  found  to  be  within  ±  1%  in  nearly  all 
cases. 


3  GEOMETRICALLY  CONSTRAINED  DESIGNS 

In  addition  to  the  baseline  design,  three  additional  design  candidates  are  developed  using  for  the 
optimization  problem  as  described  in  the  proceeding  sections.  All  benchmark  candidate 
optimization  problems  are  solved  using  an  active  set  algorithm.  The  non-linear  programming 
algorithm  is  described  by  Powell  [15]  and  incorporated  in  Matlab. 

3.1  Flat  Plate  Design 

The  baseline  for  all  candidate  design  comparisons  is  a  plate  in  which  all  nodal  locations  are  initially 
zero.  The  objective  function  is  a  single  variable  function  in  which  only  the  thickness  c  of  the  plate  is 
varied  until  the  constraints  are  satisfied.  The  optimization  problem  for  the  flat  plate  candidate 
design  is: 

find  c 

minimize  M(c)  =  pAc 

subject  to  Pc(c )  -  Pcmax  ^  0  (4) 

z  =  0 

7.5  mm  <  c  <  50  mm 

The  thickness  of  the  plate  is  computed  as  constant  at  every  node.  Initial  z-coordinate  values  are 
uniform  for  the  initial  surface,  i.e.,  all  pre-blast  nodal  locations  z  —  0.  The  Pcmax  value  allows  the 
plate  to  deflect  a  small  amount  while  still  satisfying  the  constraints,  but  the  optimized  design 
exhibits  a  relatively  large  increase  in  mass  in  order  to  absorb  the  blast  wave  while  ensuring  that  the 
deflection  constraint  is  satisfied.  This  result  illustrates  the  problem  associated  with  the  up- 
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armoring  of  vehicles  using  only  thick,  flat  plate  structures:  the  large  increase  in  vehicle  mass 
associated  with  ensuring  an  appropriate  level  of  protection  is  often  enough  to  cause  other  vehicle- 
performance  issues.  The  mass  characteristics  of  all  proceeding  designs  are  compared  to  this  result 
as  a  measure  of  design  performance.  The  results  of  the  active  set  method  optimization  are  thickness 
c  =  29.2  mm  with  a  corresponding  mass  M  =  226.4  kg. 


3.2  Pyramid  Profile  Design 


There  are  well-documented  studies  that  demonstrate  the  ability  of  v-shaped  plate  structures  to 
deflect/absorb  the  energy  from  a  blast  event.  This  candidate  design  is  meant  to  emulate  that  V- 
shaped  profile  in  rendering  a  symmetric  pyramidal  design  for  optimization.  As  the  base  of  the 
pyramid  is  fixed  at  all  constrained  nodes,  the  analytical  description  for  the  plate  structure  is  a 
function  of  the  height  of  the  structure.  The  number  of  design  variables  for  this  method  is  increased 
to  two:  the  height  of  the  structure  and  the  thickness  of  the  plate  are  optimized.  After  the  profile 
generation  procedure,  the  plate  takes  on  a  pyramid  profile  where  the  height  h  and  thickness  c  of 
the  plate  are  allowed  to  vary  and  is  optimized  with  the  objective  function  as  a  two  variable  function 
as 


find  c,  h 

minimize  M(c,  h)  =  2  pc  Jh  +  I/4 

subject  to  Pc(c,/l)  -  Pcmax  ^  0 
h 

z  -  — —  x  +  h  =  0 
0.5 

7.5  mm  <  c  <  50  mm 
0  mm  <  h  <  120  mm 


(5) 


Preliminary  results  showed  that  the  pyramid  tends  toward  the  maximum  allowable  height; 
therefore,  the  maximum  height  of  the  design  was  chosen  as  120  mm  for  the  optimization  in  order  to 
ensure  that  the  envelope  constraints  were  not  violated.  Convergent  results  demonstrate  an 
increase  in  height  to  the  boundary  of  the  design  domain  while  exhibiting  a  thickness  less  than  half 
of  that  of  the  baseline  design,  resulting  in  a  structure  of  significantly  smaller  mass. 
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Figure  6:  Finite  Element  model  showing  the  optimized  pyramid  profile  design  before  (top)  and  after  a 
blast  event  (bottom).  Convergent  numerical  results:  thickness  =  12.4  mm,  height  =  120  mm,  and  mass 

=  99.1  kg. 

As  can  be  seen  Figure  6,  the  blast  deforms  the  structure  in  such  a  way  that  the  pyramid  structure  is 
nearly  flattened,  but  the  deflection  of  the  plate  in  the  z-direction  is  within  the  allowable  tolerance. 

3.3  Gaussian  Function  Design 


Continued  investigation  of  designs  with  more  complex  curvature  led  to  the  development  of  a  plate 
structure  with  topography  that  conforms  to  the  surface  described  by  a  Gaussian  function  of  two 
variables.  The  common  single  variable  equation  for  a  Gaussian  function  is  given  as 

C x-b )21 


f(x)  =  a  exp 


2  c2 


(6) 


where  a ,  b,  and  c  are  all  constants.  In  three-dimensional  space,  given  a  mesh  grid  of  x  and  y  values, 
we  can  define  a  surface  that  is  described  by  the  equation 


f(x,y)  =  S  exp 


(x  -  x0)2  |  (y  -  y0)2 


(7) 


where  S  is  the  height  of  the  surface,  a  is  a  constant  related  to  the  spread  of  the  bulge  in  the  center. 
The  curve  is  centered  as  point  (x0,  y0).  This  allows  us  to  optimize  a  plate  structure  that  can  take  on 
complex  curvatures,  while  allowing  the  use  of  the  common  optimization  problem  used  throughout 
and  only  increasing  the  dimensional  order  of  the  design  problem  by  one. 


After  the  profile  generation  procedure,  the  plate  takes  on  a  Gaussian  profile  where  a  and  8  and 
thickness  of  the  plate  are  allowed  to  vary  and  is  the  design  is  optimized  with  the  objective  function 
as  a  function  of  the  three  variables  c,  a,  and  8  to  as 
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find 

minimize 
subject  to 


c,  a,  8 

M(c,  a,  8)  =  p  c  A(  a,  8) 
Pc(c,  a,  8)  - 

Pc  max  —  0 


z  —  8  exp 


=  0 


7.5  mm  <  c  <  50  mm 
0  <  a  <  1 
0  <  8  <  0.12 


(8) 


The  optimized  shape  before  and  after  blast  events  is  shown  in  Figure  7,  The  design  exhibits 
deformational  characteristics  similar  to  the  pyramid  profile  design  in  that  is  nearly  flattened  by  the 
blast  pressure  while  deflecting/absorbing  enough  of  the  energy  to  allow  only  a  small  amount  of 
deflection  in  the  negative  z  —direction. 


Figure  7;  FE  model  showing  the  optimized  Gaussian  profile  design  before  (top)  and  after  a  blast 
event  (bottom).  Convergent  numerical  results:  thickness  =  11.8  mm ,  a  =  0.099,  8  =  0.12,  mass  =  94.3 

kg. 

The  mass  of  the  convergent  design  is  significantly  lower  than  that  of  the  baseline  design  ( -58.4  %) 
and  improves  upon  the  mass  reduction  exhibited  by  the  pyramid  profile  design  by  a  small  margin 
( -4.84  %).  In  this  case,  an  increase  in  the  number  of  design  variables  resulting  in  a  more  complex 
shape  also  yielded  a  more  successful  design. 


3.4  Polynomial  Function  Design 


In  the  interest  of  creating  more  complex  curved  surfaces  that  could  be  controlled  through  the 
variation  of  only  a  few  variables  and  increase  the  dimension  of  the  problem  further,  a  method  for 
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generating  plate  structures  that  take  on  a  3rd  order  polynomial  curve  was  developed.  Given  the 
design  space  as  stated  above,  a  3rd  order  polynomial  curve  is  assumed  to  begin  at  some  height  (h)  at 
the  center  of  the  plate  and  end  at  a  height  of  zero  at  some  distance  (R)  from  the  center.  The 
polynomial  function  can  then  be  used  to  generate  a  surface  by  rotating  the  curve  about  the  z-axis.  It 
should  be  noted  that  this  method  does  not  make  use  of  the  entire  square  design  domain  of  the  plate 
but  instead  utilizes  a  radial  design  domain  extending  from  the  center  of  the  plate. 

Staring  with  a  third  order  polynomial  equation 

fix)  =  C0  +  Ctx  +  C2x2  +  C3x3  (9) 

we  apply  our  boundary  conditions  stated  above  as  equality  constraints  in  the  optimization  problem, 
which  results  in  a  plate  design  method  as  a  function  of  five  variables:  the  four  polynomial 
coefficients  and  the  thickness  of  the  plate.  This  in  turn  increases  the  dimension  of  the  problem  to 
five  and  the  plate  takes  on  a  polynomial  profile,  which  can  be  optimized  as  a  function  of  five 
variables: 

find  c,  Cq,  C2,  C 3 

minimize  M(c,  C0,  Ct,  C2,  C3)  =  p  A(C0,  C1,  C2,  C3)  c 
subject  to  Pc(c,  C0,  C1,  C2,  C3)  -  Pcmax  ^  0 

x?  +  yf  =  rf 

z-  C0+  Ctr  +  C2r2  +  C3r3  =  0  (10) 

h  =  120  mm 
7.5  mm  <  c  <  50  mm 
C0  h  <  0 

C0  +  RCt  +  R2C2  +  R3C3  =  0 

The  optimized  shape  before  and  after  blast  events  is  shown  in  Figure  8  and  similarly  to  the  previous 
candidate  designs,  the  plate  structure  exhibits  deformation  characteristics  in  which  the  structure  is 
nearly  flattened  by  the  blast,  while  allowing  only  a  small  amount  of  deflection  in  the  z-direction 
below  the  datum  plane,  within  the  tolerance  of  the  penetration  constraint.  The  optimized  design  is 
an  improvement  over  the  baseline  design  in  terms  of  mass  reduction,  but  demonstrates  an  increase 
in  comparison  to  the  pyramid  profile  and  Gaussian  profile  designs  due  to  an  increase  in  thickness  in 
the  converged  results.  It  is  thought  that  this  increase  in  thickness  is  in  compensation  for  a  feature  of 
this  design  nit  present  in  the  previous  two  design  methods  mentioned  here.  Due  to  the  fact  that  the 
design  is  generated  as  a  surface  of  revolution,  which  is  inherently  circular,  while  the  shape  of  the 
plate  as  described  by  the  domain  is  square;  the  whole  design  domain  is  not  used  on  the  generation 
of  the  plate  profile.  In  turn,  this  creates  an  inherent  weakness  in  the  plate  structure  in  the  portion  of 
the  plate  beyond  the  area  described  by  the  surface  of  revolution  -  i.e.  the  flat  portions  of  the  plate 
outside  of  the  bulging  center  portion. 
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Figure  8:  Finite  Element  model  showing  the  optimized  polynomial  function  surface  design  before 
(top)  and  after  a  blast  event  (bottom).  Convergent  numerical  results:  thickness  =  14.1  mm,  C0  = 
0.112,  C0  =  -0.332,  C0  =  0.642,  C0  =  -0.859,  mass  =  113.1  kg. 

3.5  Trigonometric  Function  Design 

The  trigonometric  function  design  was  developed  in  an  effort  to  explore  a  design  whose  surface 
exhibits  a  more  complex  curvature  than  a  3rd  order  polynomial  without  substantially  raising  the 
dimension  of  the  design  problem.  Although  the  polynomial  design  did  not  demonstrate  an  increase 
in  performance  correlated  to  an  increase  in  problem  dimension,  this  design  allows  a  dramatic 
increase  in  curvature  complexity  with  only  a  marginal  increase  in  the  problem  dimension.  While 
the  shape  generation  process  is  very  similar  to  the  one  used  in  the  polynomial  function  design,  the 
surface  of  revolution  is  based  on  the  sum  of  a  set  of  oscillating  functions,  namely  cos  and  sin. 
Through  the  definition  of  a  few  coefficients,  this  function  allows  for  curvature  that  could  only  be 
represented  by  a  very  high-order  polynomial  function.  Similar  to  the  polynomial  design  method 
described  above,  this  method  does  not  make  use  of  the  entire  square  design  domain  of  the  plate  but 
instead  utilizes  a  radial  design  domain  extending  from  the  center  of  the  plate. 

Given  the  design  space  as  stated  above,  a  curve  is  assumed  to  begin  at  some  height  (h)  at  the  center 
of  the  plate  and  end  at  a  height  of  zero  at  some  distance  (R)  from  the  center.  We  define  /(x)  as  the 
sum  of  a  cos  and  sin  function  with  coefficients  to  allow  customization  of  the  curve  within  the  design 
domain: 

/(x)  =  CqCOsClX  +  C2s\nC3x  (11) 

After  some  development,  the  coefficients  are  chosen  to  allow  the  designer  to  control  overall  height 
and  width  of  the  structure  and  allow  for  some  amplitude  and  frequency  modulation  of  the  second 
term  to  generate  a  wave  shape.  This  would  allow  for  the  creation  of  plate  structures  of  much 
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greater  complexity  than  a  3rd  order  polynomial.  The  coefficients  in  the  equation  above  become  the 
design  variables  and  the  trigonometric  function  for  design  generation  is: 


fix)  =  %  cosnx  + 


ia 21  a22) 


R 


x  +  a2 1 


sin  ^ 


(A  ~  A) 

R 


x  + 


f 2) 


2nx 


(12) 


This  method  exhibits  an  increase  in  problem  dimension  to  six  and  the  plate  takes  on  a 
trigonometric  function  profile,  which  can  be  optimized  via  gradient-based  methods  as  a  function  of 
six  variables.  The  optimization  problem  is  stated  as 

find  c,  a1(  a2i,  a22>  fi>  f2 

minimize  Mic,a1,a21,a22,f1,f2)  =  P  A(c,  a1,  a21,  a22,A,A)  c 


subject  to 


Pdc,  al>a21>a22>fl’f2)  ~  Pc  max  —  0 


xf  +  yf  =  rf 


z  —  at  cos  nr  + 


ia 21  ~  a22) 


R 


r+  a21 


sin 


(/2  -  A) 


R 


r  +  f2  1  2nr 


=  0 


7.5  mm  <  c  <  50  mm 
h  =  120  mm 
0  <  ax  <  h—  a21 
0  mm  <  a21  <  20  mm 
0  mm  <  a22  <  5  mm 
0  <  A  <  50 
0  <  /2  <  50 


(13) 


3.5.1  Mesh  refinement  and  two-stage  optimization 


The  complex  curvature  of  the  surfaces  generated  by  the  trigonometric  function  design  were  not 
able  to  be  rendered  with  accuracy  in  the  standard  26  x  26  plate  used  throughout  this  study.  The 
design  can  exhibit  many  small  rippling  curvatures,  which  cannot  be  captured  by  a  course  mesh.  As 
shown  in  Figure  9,  an  increase  in  the  level  of  mesh  refinement  is  necessary  to  obtain  valid 
experimental  results. 


Figure  9:  Detail  of  plate  showing  mesh  refinement  to  capture  detailed  curvature  of  trigonometric 
function  design  method.  100x100  element  mesh  refinement  (left)  and  200x200  element  mesh 

refinement  (right)  are  shown. 


Through  a  mesh  convergence  study  it  was  found  that  a  200  x  200  mesh  discretization  of  the  design 
domain  was  adequate  to  achieve  accurate  results.  This  mesh  refinement  comes  at  the  cost  of 
dramatically  increased  computational  time  however,  as  the  number  of  elements  increases  from  625 
to  40000.  This  increase  correlates  to  an  increase  in  FEA  simulation  time  from  approximately  5 
seconds  to  over  10  minutes  in  some  cases.  When  the  bounds  on  the  design  variables  are  large,  the 
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number  of  function  calls  necessary  to  find  the  gradient  for  optimization  can  be  high,  resulting  in  a 
convergence  time  for  this  optimization  problem  that  is  unreasonably  high.  In  order  to  overcome 
this  problem,  the  optimization  is  carried  out  in  two  stages  of  fidelity.  A  model  in  which  the  design 
domain  is  discretized  into  100  x  100  elements,  considered  a  low-fidelity  model  with  a  relatively 
large  error,  is  used  to  narrow  the  bounds  to  within  reasonable  spans.  After  the  span  of  bounds  has 
been  reduced,  a  high-fidelity  model  is  used  to  find  the  convergent  design.  This  process  resulted  in 
an  accurate,  convergent  design  at  a  fraction  of  the  computational  cost. 

The  convergent  design,  shown  in  Figure  10,  exhibits  some  amplitude  and  frequency  modulation  of 
the  second  trigonometric  term,  and  the  plate  takes  on  the  appearance  of  raised  concentric  ripples 
emanating  from  the  center  of  the  structure. 


Figure  10:  Finite  Element  model  showing  the  optimized  trigonometric  function  surface  design  before 
(top)  and  after  a  blast  event  (bottom).  Convergent  numerical  results:  thickness  =  8.31  mm ,  at  = 
0.1045 ,  a21  =  0.01549 ,  a22  =  0.0000,  f1  =  44.794,  f2  =  34.000,  mass  =  81.2  kg. 


This  design  method  yields  the  greatest  reduction  in  mass  over  baseline  (-64.1  %)  and  a  substantial 
reduction  in  mass  over  the  Gaussian  design,  the  next  best  performing  method  (-14.0  %).  Despite  the 
increase  in  the  complexity  of  the  curvature,  the  reduction  in  mass  is  due  to  a  thinner  structure  than 
was  found  in  previous  results.  In  this  case,  an  increase  in  the  order  of  the  optimization  problem 
allowed  for  a  more  pronounced  relaxation  of  the  design,  and  resulted  in  a  higher  performance, 
more  complex  structure. 
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4 


FREE  SHAPE  DESIGNS 


4.1  Inverted  Profile  Design 

There  has  been  some  interest  in  creating  a  plate  that  takes  on  an  inversion  of  the  topography  of  a 
flat  plate  after  a  blast  event.  The  reasoning  being  that  the  complex  blast  pressure  distribution 
would  then  determine  the  shape  of  the  plate  and  increase  the  z-coordinate  at  each  proportional  to 
the  magnitude  of  the  loading.  The  basic  procedure  used  to  generate  the  topography  for  this  plate  is 
as  follows: 

•  Expose  a  flat  plate  of  some  thickness  c  to  a  blast  event  of  time  t  —  tc  in  order  to  capture  the 
maximum  deflection  caused  by  the  blast. 

•  Invert  the  shape  about  the  xy-plane  for  all  nodal  locations. 

•  Scale  all  z  values  to  remain  within  the  design  space  (if  necessary). 

Through  the  execution  of  this  procedure,  the  topography  of  the  plate  is  determined  without  the 
application  of  any  rigorous  analytical  description.  This  inversion  procedure  allows  the  designer  to 
create  complex  shapes  proportional  to  the  loading  condition  without  the  need  for  complicated 
analytical  descriptions  that  would  not  be  feasibly  solvable  by  traditional  gradient-based 
optimization  methods.  After  the  inversion  procedure,  the  topography  of  the  plate  is  fully  specified, 
and  the  design  is  optimized  via  the  active-set  algorithm  used  to  find  the  convergent  results  for  the 
geometrically  constrained  designs  in  the  previous  sections.  The  objective  function  is  a  single 
variable  function  where  only  the  thickness  of  the  plate  structure  is  allowed  to  vary: 

find  c 

minimize  M(c)  =  pAc  „ 

(14) 

subject  to  Pc(c)  -  Pcmax  ^  0 
7.5  mm  <  c  <  50  mm 

The  optimized  design,  shown  in  Figure  11,  exhibits  a  performance  increase  over  the  baseline  design 
in  terms  of  mass  reduction  very  similar  to  the  geometrically  constrained  designs  but  with  an 
optimization  problem  of  order  one,  whose  solution  can  be  found  with  much  less  computation.  In 
addition,  it  utilizes  the  entire  square  area  of  the  design  domain  and  could  easily  be  applied  to 
irregularly  shaped  design  domains  with  little  difficulty. 
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Figure  11:  Finite  Element  model  showing  the  optimized  inverted  profile  design  before  (top)  and  after 
a  blast  event  (bottom).  Convergent  numerical  results:  thickness  =  13.6  mm ,  mass  =  105.8  kg. 


4.2  HCA  Topography  Design 

The  HCA  topography  design  is  an  alternative  free  shape  design  method,  which  allows  for  the 
generation  of  complex  curved  structures  without  the  need  for  complex  analytical  description.  The 
HCA  algorithm  considers  the  nodal  z-coordinates  of  all  unconstrained  nodes  in  the  design  domain 
at  every  iteration.  Thus,  the  objective  function  is  of  dimension  N,  where  N  is  the  total  number  of 
free  nodes  in  the  design  domain.  The  convergence  behavior  is  such  that  the  shape  depth  is 
increased  incrementally  at  each  iteration  of  the  control-based  HCA  algorithm  until  the  penetration 
constraint  is  satisfied.  The  thickness  of  the  plate  is  fixed,  and  the  objective  function  M,  is  purely  a 
function  of  the  plate  geometry.  The  optimizer  does  not  evaluate  the  mass  at  every  iteration,  but 
simply  iterates  until  the  constraints  are  met.  Given  an  initially  flat  surface  (all  nodal  locations  z  = 

0)  and  a  fixed  elemental  thickness  (c  =  q),  the  HCA  algorithm  is  applied  to  the  design  domain 
iteratively  until  convergence  criteria  as  fixed  by  the  designer  are  met. 

After  preliminary  analysis,  a  thickness  of  q  =  12.2  mm  was  chosen  in  order  to  keep  the  structure  as 
lightweight  as  possible  while  utilizing  the  entire  design  domain  without  violating  deflection 
constraints.  The  final  design  obtained  from  the  nodal  HCA  algorithm  is  given  in  Figure  12  below, 
with  the  resulting  parameter  information  given. 
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Figure  12:  Finite  Element  model  showing  the  optimized  HCA  profile  design  before  (top)  and  after  a 
blast  event  (bottom).  Numerical  results:  thickness  =  12.2  mm ,  structure  height  =11.9  mm,  mass  = 

97.8  kg. 


To  demonstrate  the  dependence  of  the  converged  solution  on  user  inputs,  a  second  design  is 
generated  using  an  initial  starting  thickness  of  q  =  0.015  m.  This  second  structure  satisfies  the 
same  penetration  constraint  with  a  lower  height,  while  exhibiting  an  increase  in  mass.  In  other 
words,  the  plate  is  designed  to  be  thicker  and  have  more  mass  with  the  trade-off  that  it  takes  up 
less  space  in  the  vertical  direction,  demonstrating  the  high-level  of  control  the  designer  can 
maintain  in  the  customization  of  the  structure.  Table  2  shows  the  comparative  results. 

Table  2:  Numerical  comparison  of  HCA  profile  designs  1  and  2. 


Trial 

Thickness 

(mm) 

Height  of 
Structure  (mm) 

Mass  (kg) 

HCA  Design  1 

12.2 

11.9 

97.8 

HCA  Design  2 

15.0 

8.61 

118.7 

While  the  geometrically  constrained  candidate  designs  can  be  optimized  with  gradient-based 
methods,  this  approach  is  not  feasible  for  the  HCA  objective  function  due  to  the  large  dimensional 
order  of  the  optimization  problem.  Convergence  criteria  must  instead  be  chosen  by  the  designer 
based  on  specific  performance  goals  and  the  application  of  the  structure.  In  addition,  the  design 
rules  are  extremely  problem  dependent  and  are  determined  through  algorithmic  testing.  Despite 
the  problems  of  demonstrating  convergence  and  the  heuristic  nature  of  the  design,  the  HCA  method 
exhibits  some  key  advantages  over  geometrically  constrained  designs:  computational  efficiency  and 
a  broad  potential  for  application  to  irregular  design  domains. 
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4.2.1  Design  algorithm 


The  hybrid  cellular  automaton  (HCA)  method  combines  the  basis  of  the  cellular  automaton  (CA) 
paradigm;  introduced  by  Stanislaw  Ulam  and  John  von  Neumann  in  1950s  [16],  and  the  theory  of 
finite  element-based  structural  optimization;  introduced  by  Lucien  Schmit  in  1960s  [17].  The  HCA 
method  presented  by  Tovar,  et  al.  [18]  incorporates  local  updating  schemes  such  as  control  rules 
(i.e.,  on-off,  proportional,  integral  and  derivative  controllers  and  ratio  techniques).  These  local  rules 
drive  a  defined  field  variable  to  an  optimum  state  or  set  point.  The  expression  for  the  field  variable 
and  the  value  of  the  set  point  are  derived  from  the  optimality  conditions  of  the  structural  design 
problem  [19].  A  proof  of  the  global  convergence  of  the  HCA  technique,  under  certain  circumstances, 
to  an  optimal  design  has  been  derived  [1]. 


Figure  13:  HCA  algorithm  process  flow  diagram  showing  plate  structure  application. 


The  HCA  algorithm  as  applied  to  the  plate  topography  problem  is  illustrated  in  Figure  13  above. 
This  methodology  assumes  that  cellular  automata  (CAs)  form  a  structure  or  design  domain,  and 
sensors  and  actuators  within  the  CAs  activate  local  formation  and  resorption  of  material.  With  a 
proper  control  strategy,  this  process  drives  the  overall  structure  to  an  optimal  topography  by 
updating  the  nodal  locations.  Using  distributed  controlled  rules,  the  optimization  problem  can  be 
stated  as 
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(15) 


find  z 

minimize  \q*  -  qmin\ 

Subject  to  l^minl  —  Pc  max  —  0 
7^  <  2  <  z U 

where  z  is  the  set  of  all  nodal  coordinate  values  associated  the  plate  topography.  In  (15)  q*  is  the 
target  deflection  threshold  value  to  be  achieved  by  every  element  in  the  structure.  The  iterative 
approach  may  be  achieved  using  a  control-based  algorithm  or  a  ratio  approach.  A  control-based 
approach  is  implemented  to  perform  topography  optimization  as  shown  in  Figure  14.  Multiple 
gains  [Ka  and  Kp]  allow  the  designer  control  over  the  rate  at  which  the  topography  of  the  structure 
is  updated. 


<r 


z 

Figure  14:  Control-based  strategy  for  updating  nodal  locations  to  optimize  topography. 


4.2.2  Computational  efficiency  and  convergence  considerations 

The  computational  efficiency  of  the  HCA  topography  design  is  rooted  in  two  key  characteristics  of 
the  design  method.  The  first  is  the  relative  insensitivity  to  mesh  refinement.  Due  to  the  nature  of 
the  shape  produced  by  HCA  topography  a  fairly  coarse  mesh  can  be  used  to  generate  structures,  as 
opposed  to  some  of  the  more  complex  geometrically  constrained  design  methodologies,  which 
require  a  higher  level  of  mesh  refinement  for  modeling  accuracy.  A  mesh  convergence  study 
performed  on  the  HCA  method  demonstrated  nearly  identical  results  for  models  containing  625 
elements  and  40000  elements  in  terms  of  mass  and  deflection. 


UNCLASSIFIED  Dist  A.  Approved  for  public  Release. 


22 


Table  3:  HCA  mesh  convergence  study  results. 


Grid  size 

Number 

of  elements 

Mass  of 
converged 
design  (kg) 

25x25 

625 

97.8 

100x100 

10000 

98.1 

200x200 

40000 

98.1 

In  addition  to  the  relatively  low  mesh  requirements,  the  HCA  method  is  in  general  more 
computationally  efficient  in  that  it  only  requires  one  function  call  per  iteration,  independent  of  the 
number  of  design  variables.  Gradient-based  optimization  methods  can  require  many  function  calls, 
in  some  cases  hundreds  per  iteration  to  determine  the  gradient  and  hessian  matrices  for  higher- 
dimension  optimization  problems.  In  the  case  of  structure  design,  these  multiple  function  calls  to 
FEA  simulation  can  be  very  computationally  expensive.  This  benefit  is  augmented  by  the  controls- 
based  optimization  routine,  which  allows  the  designer  to  control  the  rate  of  convergence  through 
gain  adjustments  in  the  control  rule. 

4.2.3  Application  to  irregular  design  domains 

The  other  primary  advantage  of  the  HCA  topography  has  over  the  geometrically  constrained 
designs  is  its  potential  for  much  broader  application  to  irregular  design  domains.  While  some  of  the 
geometrically  constrained  designs  are  dependent  on  symmetric  bounds  (e.g.  the  pyramid  profile 
design)  and  others  fail  to  take  advantage  of  the  entire  design  domain  due  to  their  inherently 
circular  footprint  (e.g.  the  polynomial  and  trigonometric  function  designs),  the  HCA  can  be  applied 
without  difficulty  to  a  wide  range  of  design  domains  under  vastly  different  loading  and  support 
conditions.  Figure  15  below  demonstrates  the  adaptation  of  the  HCA  method  to  an  irregular  design 
space  meant  to  simulate  the  under-armoring  of  a  vehicle. 


Figure  15:  HCA  topography  design  in  application  to  under-vehicle  design  domain.  The  domain  is 
shown  from  the  bottom  of  the  vehicle  (left)  and  the  convergent  design  is  shown  in  application  from 

the  right  (right). 
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5  COMPARATIVE  NUMERICAL  RESULTS 


The  numerical  results  of  the  convergent  designs  are  divided  into  two  tables  in  order  to  compare 
only  designs  operating  within  the  same  design  domain.  Table  4  contains  performance  values  for  the 
radial  design  domain  methods.  Table  5  contains  performance  values  for  the  full  design  domain 
methods.  In  both  design  domain  cases,  radial  and  full,  there  is  a  trend  of  increased  performance  in 
terms  of  mass  reduction  as  the  number  of  design  variables  increases. 


Table  4:  Comparative  numerical  results:  radial  design  domain 


Method 

No.  of 
design 
variables 

No.  of 

iterations 

No.  of 

function 

calls 

Mass  of  converged 
design  (kg) 

Mass  reduction 
from  baseline  (%) 

Flat  Plate  Design 

1 

18 

71 

226.4 

0.0 

Polynomial  Function  Design 

5 

173 

1398 

113.1 

-50.0 

Trigonometric  Function  Design 

6 

272* 

1723* 

81.2 

-64.1 

*  These  totals  are  a  sum  of  the  iterations  necessary  for  the  high  and  low  fidelity  models  to  converge 

Table  5:  Comparative  numerical  results:  full  design  domain 

Method 

No.  of 
design 
variables 

No.  of 

iterations 

No.  of 

function 

calls 

Mass  of  converged 
design  (kg) 

Mass  reduction 
from  baseline  (%) 

Flat  Plate  Design 

1 

18 

71 

226.4 

0.0 

Pyramid  Profile  Design 

2 

17 

121 

99.1 

-56.2 

Gaussian  Function  Design 

3 

9 

82 

94.3 

-58.4 

Inverted  Profile  Design 

n/a** 

n/a 

n/a 

105.8 

-53.3 

HCA  Topography 

112 

112 

97.8 

-56.8 

**  The  procedure  of  inverting  the  shape  does  not  apply  to  the  computational  efficiency  of  the  design  method 
***  N  is  equal  to  the  number  of  free  nodes  in  the  design  domain 


As  was  stated  previously,  the  validation  of  HCA  concept  designs  can  be  extremely  problematic,  but 
an  attempt  is  made  here  to  accomplished  validation  through  direct  comparison  to  benchmark 
results.  The  summary  of  numerical  results  as  an  improvement  of  the  baseline  design  is  given  in 
Table  5.  In  every  case,  the  penetration  for  the  converged  design  is  zero.  The  comparative  numerical 
results  for  the  HCA  method  demonstrate  a  substantial  performance  improvement  over  the  baseline 
design  in  terms  of  mass  reduction  and  performance  approximately  equal  to  or  better  than  all 
converged  candidate  designs,  with  the  exception  of  the  trigonometric  function  design  method. 
These  results  demonstrate  the  success  of  the  HCA  topography  design  in  generating  structures  of 
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equal  and  in  some  cases  superior  to  geometrically  constrained  designs  while  using  fewer 
computational  resources  and  allowing  the  designer  a  higher  degree  of  control. 


6  CONCLUSIONS 

In  this  work,  the  topography  optimization  of  isotropic  steel  plate  structures  to  mitigate  blast 
loading  is  accomplished.  Several  design  strategies  are  examined,  falling  under  two  categories: 
geometrically  constrained  designs,  whose  shape  is  bound  to  analytical  definition,  and  free  space 
designs,  whose  shape  is  allowed  to  change  freely  under  the  control  of  the  optimization  algorithm.  In 
the  optimization  problem  applied  to  all  designs,  the  mass  of  the  structure  is  minimized  while 
deflection  constraints  are  met.  Through  this  comparison  study,  it  is  possible  to  draw  some 
conclusions  about  the  success  of  the  various  methods  to  produce  good  candidates  for  blast 
protection  design. 

Through  optimization  of  topography,  all  designs  exhibit  an  increase  in  performance  in  terms  of 
mass  reduction  over  baseline  while  maintaining  an  allowable  amount  of  deflection.  Optimum 
designs  found  through  an  expansion  of  the  order  of  the  optimization  problem  from  one  to  six 
demonstrated  a  general  trend  of  increased  performance  as  the  domain  constraints  are  relaxed  and 
the  design  is  able  to  assume  more  complex  analytical  shapes.  Numerical  results  of  convergent 
designs  based  on  trigonometric  functions  demonstrate  a  reduction  in  mass  of  64%  over  baseline. 
This  result  exceeds  all  other  performance  results  for  the  various  designs  and  makes  the 
trigonometric  function  design  a  good  candidate  for  further  study,  especially  the  potential  to  apply 
the  methodology  to  irregular  design  spaces  and  the  sensitivity  of  the  design  to  variations  in  loading 
conditions.  Future  research  is  indicated  in  the  development  of  a  robust  design  rooted  in  this 
trigonometric  function  design  method. 

It  is  problematic  to  prove  a  global  optimum  with  HCA  and  other  free  space  methods,  as  their 
method  of  optimization  is  not  based  on  KKT  conditions.  However,  designs  generated  via  the  Hybrid 
Cellular  Automata  (HCA)  topography  method  resulted  in  structures  whose  performance  paralleled 
or  exceeded  most  geometrically  constrained  designs  despite  the  fact  that  they  can  be  created  with 
fewer  computational  resources  and  adapted  easily  to  irregular  design  domains.  Under  well-chosen 
control  gains,  the  HCA  method  can  produce  designs  similar  in  quality  to  that  of  designs  produced  by 
gradient-based  optimization  methods  using  hundreds  of  fewer  function  calls. 

The  primary  disadvantage  of  this  type  of  empirical  model  for  blast  loading  is  the  absence  of 
reflective  waves  from  surrounding  surfaces,  i.e.,  the  ground,  buildings,  etc.  Future  studies  should 
include  MM-ALE  or  similar  solver  to  account  for  fluid-material  interactions  and  include  those 
effects.  Much  attention  has  been  paid  recently  to  the  development  of  sandwich  structures  as  a 
multi-material  approach  to  blast  mitigation  problems  and  some  have  demonstrated  significant 
reduction  of  stress  amplitude,  especially  through  the  use  of  honeycomb  structures.  Future  work 
will  incorporate  concurrent  material  and  shape  optimization,  and  if  possible  will  include  the 
application  of  the  design  candidates  examined  here  to  multiple-material  design. 


UNCLASSIFIED  Dist  A.  Approved  for  public  Release. 


25 


7  ACKNOWLEDGEMENTS 


This  material  is  based  upon  work  supported  by  the  U.S.  Army  TACOM  Life  Cycle  Command  under 
Contract  No.  W56HZV-08-C-0236,  through  a  subcontract  with  Mississippi  State  University  and  the 
University  of  Notre  Dame,  and  was  performed  for  the  Simulation  Based  Reliability  and  Safety 
(SimBRS)  research  program.  Any  opinions,  findings,  conclusions,  and  recommendations  expressed 
in  this  paper  are  those  of  the  writers  and  do  not  necessarily  reflect  the  views  of  the  sponsors. 


8  REFERENCES 

1.  Goetz,  J.,  et  al.,  Two-material  optimization  of  plate  armour  for  blast  mitigation  using 
hybrid  cellular  automata.  Engineering  Optimization,  2012.  44(8]:  p.  985-1005. 

2.  Argod,  V.,  et  al.,  Shape  optimization  of  solid  isotropic  plates  to  mitigate  the  effects  of 
air  blast  loading.  Mechanics  Based  Design  of  Structures  and  Machines,  2010.  38(3): 
p.  362-371. 

3.  Wang,  H.,  et  al.  Function-Oriented  design  of  innovative  composite  materials  for  high 
speed  impact,  in  Materials  Science  and  Technology  Conference  and  Exhibition,  MS  and 
T'07  -  "Exploring  Structure,  Processing,  and  Applications  Across  Multiple  Materials 
Systems",  September  16,  2007  -  September  20,  2007.  2007.  Detroit,  MI,  United  states: 
Curran  Associates  Inc. 

4.  Jiang,  D.,  et  al.,  Innovative  composite  structure  design  for  blast  protection,  in  SAE 
World  Congress2007,  US  Army  RDECOM-TARDEC:  6501  E  11  Mile  Rd  Warren  MI 
48397-5000. 

5.  Tan,  P.J.,  et  al.,  Dynamic  compressive  strength  properties  of  aluminium  foams.  Part  II  - 
'shock'  theory  and  comparison  with  experimental  data  and  numerical  models.  Journal 
of  the  Mechanics  and  Physics  of  Solids,  2005.  53(10):  p.  2206-2230. 

6.  Hanssen,  A.G.,  et  al.  Close-range  blast  loading  of  aluminium  foam  panels:  A  numerical 
study,  in  IUTAM  Symposium  on  Mechanical  Properties  of  Cellular  Materials,  September 
1 7, 2007  -  September  20,  2007.  2009.  Cachan,  France:  Springer  Verlag. 

7.  Randers-Pehrson,  G.  and  K.A.  Bannister,  Airblast  loading  model  for  DYNA2D  and 
DYNA3D,  1997:  Army  Research  Laboratory. 

8.  Kuhl,  A.L.,  Dynamics  of  detonations  and  explosions  -  explosion  phenomena,  1991, 
American  Institute  of  Aeronautics  and  Astronautics:  Washington,  D.C. 

9.  Hirt,  C.W.,  A.A.  Amsden,  and  J.L.  Cook,  An  arbitrary  lagrangian-eulerian  computing 
method  for  all  flow  speeds.  Journal  of  Computational  Physics,  1974. 14:  p.  227-253. 

10.  Joynt,  V.P.,  Mine  Resistant  Armored  Vehicle,  2008,  Force  Protection  Industries,  Inc. 
(Ladson,  SC,  US):  United  Sates. 

11.  Tan,  H.,  Vehicle  Hull  Shape  Optimization  for  Minimum  Deformation  Under  Blast 
Loading  2010. 

12.  Troyano,  L.F.,  Bridge  Engineering:  A  Global  Perspective.  2003:  Thomas  Teflord. 

13.  Osserman,  R.,  How  the  Gateway  Arch  Got  its  Shape.  NEXUS  NETWORK  JOURNAL, 
2010. 12(2):  p.167-189. 

14.  Ngo,  T.,  et  al.,  Blast  loading  and  blast  effects  on  structures  -  An  overview.  Electronic 
Journal  of  Structural  Engineering,  2007.  7:  p.  76-91. 


UNCLASSIFIED  Dist  A.  Approved  for  public  Release. 


26 


15.  Powell,  M.J.D.,  A  Fast  Algorithm  for  Nonlinearly  Constrained  Optimization 
Calculations.  Numerical  Analysis,  1978.  Vol.  630. 

16.  Wolfram,  S.,  A  New  Kind  of  Science 2002,  Champaign,  IL:  Wolfram  Media.  1197. 

17.  Schmit,  L.A.,  Some  approximation  concepts  for  structural  synthesis.  AIAA  Journal, 
1974. 12:  p.692-699. 

18.  Tovar,  A.,  et  al.,  Topology  Optimization  Using  a  Hybrid  Cellular  Automaton  Method 
with  Local  Control  Rules.  ASME  Journal  of  Mechanical  Design,  2006. 128(6):  p.  1205- 
1216. 

19.  Tovar,  A.,  et  al.,  Optimality  Conditions  of  the  Hybrid  Cellular  Automata  for  Structural 
Optimization.  AIAA  Journal,  2007. 45(3):  p.  673-683. 


UNCLASSIFIED  Dist  A.  Approved  for  public  Release. 


27 


